AD-A228  590 


NASA  Contractor  Report  187461 
ICASE  Report  90-78 


ICASE  OTIC  FILE  CORY 


THE  ROUTE  TO  CHAOS  FOR  THE 
KURAMOTO-SIVASHINSKY  EQUATION 


Demetrios  T.  Papageorgiou 
Yiorgos  S.  Smyrlis 


Contract  No.  NAS  1-18605 
October  1990 

Institute  for  Computer  Applications  in  Science  and  Engineering 
NASA  Langley  Research  Center 
Hampton,  Virginia  23665-5225 

Operated  by  the  Universities  Space  Research  Association 


NASA 

National  Aeronautics  and 
Space  Administration 

Langley  Research  Center 

I  lampion,  Virginia  23665  5225 


»Q  U.  9  012 


The  route  to  chaos  for  the  Kuramoto-Sivashinsky  equation.1 

Demetrios  T.  Papageorgiou 
Department  of  Mathematics  and 
Center  for  Applied  Mathematics  and  Statistics 
New  Jersey  Institute  of  Technology 
Newark,  New  Jersey  07102 

Yiorgos  S.  Smyrlis 
Department  of  Mathematics 
University  of  California  at  Los  Angeles 
Los  Angeles,  California  90024-1555 


ABSTRACT 

We  present  the  results  of  extensive  numerical  experiments  of  the  spatially  periodic  ini¬ 
tial  value  problem  for  the  KuramotoSivashinsky  equation.  Our  concern  is  with  the  asymp¬ 
totic  nonlinear  dynamics  as  the  dissipation  parameter  decreases  and  spatio-temporal  chaos  sets 
in.  To  this  end  the  initial  condition  is  taken  to  be  the  same  for  all  numerical  experiments  (a 
single  sine  wave  is  used)  and  the  large  lime  evolution  of  the  system  is  followed  numerically. 
Numerous  computations  were  performed  to  establish  the  existence  of  windows,  in  parameter 
space,  in  which  the  solution  has  the  following  characteristics  as  the  viscosity  is  decreased:  A 
steady  fully  modal  attractor  to  a  steady  bimodal  attractor  to  another  steady  fully  modal  attrac¬ 
tor  to  a  steady  trimodal  attract  to  a  periodic  (in  lime)  attractor,  to  another  steady  fully 
modal  attractor,  to  another  time-periodic  attractor,  to  a  steady  tetramodal  attractor,  to  another 
time-periodic  attractor  having  a  full  sequence  of  period-doublings  (in  the  parameter  space)  to 
chaos.  Numerous  solutions  are  presented  which  provide  conclusive  evidence  of  the  period- 
doubling  cascades  which  precede  chaos  for  this  infinite-dimensional  dynamical  system.  These 
results  permit  a  computation  of  the  lengths  of  subwindows  which  in  turn  provide  an  estimate 
for  their  successive  ratios  as  the  cascade  develops.  A  calculation  based  on  the  numerical 
results  is  also  presented  to  show  that  the  period  doubling  sequences  found  here  for  the 
Kuramoto-Sivashinsky  equation,  are  in  complete  agreement  with  Feigenbaum’s  universal  con¬ 
stant  of  4.669201609...  .  Some  preliminary  work  shows  several  other  windows  following  the 
first  chaotic  one  including  periodic,  chaotic  and  a  steady  octamodal  window;  however  the 
windows  shrink  significantly  in  size  to  enable  concrete  quantitative  conclusions  to  be  made. 
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1.  Introduction. 

The  present  study  is  concerned  with  the  long  time  behavior  of  the  one-dimensional  Kuramoto- 
Sivashinsky  (KS)  equation,  conveniently  normalized  as 

U,  +  UUX  +  Ua  +  VUxxxx  =  0  , 

(x,t)  e  RlxR+  ,  (1.1) 

u(x,Q)  =  u0(x)  ,  u{x+2K,t)  =  u(x,t)  . 

The  parameter  v  is  positive  and  takes  the  role  of  the  viscosity  of  the  system. 

The  long  time  behavior  of  the  KS  equation  is  characterized  by  the  negative  -  therefore  destabiliz¬ 
ing  -  second  order  diffusion,  the  positive  -  therefore  stabilizing  -  fourth  order  dissipation  and  the  non¬ 
linear  coupling  term.  As  soon  as  the  value  of  v  decreases  below  1  at  least  one  Fourier  frequency 
becomes  unstable  in  the  linearized  model.  Nevertheless  the  solution  does  not  grow  exponentially  fast 
as  the  linear  analysis  around  the  null  state  would  suggest  but  instead  stays  bounded  due  to  the  transfer 
of  energy  from  the  linearly  unstable  low  modes,  the  number  of  which  is  mod(v~1/2),  to  tire  linearly 
stable  high  modes.  The  previous  statement  has  been  observed  numerically  and  also  established  analyti¬ 
cally  (see  Constantin  et  al.  [8].) 

A  very  attractive  feature  of  (1.1),  is  that  even  though  it  appears  to  be  a  very  simple  partial 
differential  equation  it  produces  complicated  dynamics  in  both  space  and  time  (chaotic  behavior).  The 
equation  arises  in  a  variety  of  physical  problems,  such  as  in  chemical  physics  for  propagation  of  con¬ 
centration  waves  (Kuramoto  [24],  Kuramoto  and  Tsuzuki  [25], [26]),  plasma  physics  (Cohen  et  al.  [6]), 
flame  propagation  (Sivashinsky  [35]),  reaction  diffusion  combustion  dynamics  (Sivashinsky  [36]),  free 
surface  film-flows  (Benney  [4],  Sivashinsky  and  Michelson  [37],  Shlang  and  Sivashinsky  [33],  Hooper 
and  Grimshaw  [17]),  and  two-phase  flows  in  cylindrical  geometries  (Frenkel  et  al.  [13],  Papageorgiou 
et  al.  [31]). 

There  have  been  numerous  computational  studies  of  (1.1)  (Cohen  et  al.  [6],  Sivashinsky  and 
Michelson  [37],  Aimar  [1],  Manneville  [29],  Frisch,  She  and  Thual  [14],  Hyman  and  Nicolaenko  [18], 
Hyman,  Nicolaenko  and  Zaleski  [19],  Kevrekidis,  Nicolaenko  and  Scovel  [22],  Greene  and  Kim  [15] 
to  mention  a  few),  as  well  as  analytical  ones  (  Constantin  et  al.  [8],  Foias  et  al.  [12],  Babin  and  Vishik 
[3],  etc.).  The  overall  picture  that  emerges  is  that  as  the  viscosity  parameter  v  is  decreased,  the  large 
time  behavior  of  the  system  varies  from  steady-state  (for  example  cellular)  solutions  to  complicated 
time-oscillatory  (chaotic)  ones.  Such  behavior  has  led  to  the  conjecture  that  the  solution  behaves 
according  to  low  modal  dynamics  with  the  determining  number  of  Fourier  modes  being  proportional  to 
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v-1/2  (i.e.  Pomeau  and  Manneville  [32]).  Significant  analytical  results  have  provided  estimates  for  the 
dimension  and  the  nature  of  the  global  attractor  (see  [8], [12]),  but  the  number  of  determining  modes 
has  not  yet  been  estimated  mathematically  with  satisfactory  accuracy  and  thus  we  chose  initially  to 
Integrate  the  partial  differential  equation  directly  rather  than  work  with  a  truncated  system.  A  truncated 
system  of  coupled  nonlinear  ODE’s  for  the  Fourier  coefficients  was  also  integrated  and  the  number  of 
determining  coefficients  was  established  experimentally.  Since  the  dynamic  behavior  of  the  system  is 
of  a  low  modal  structure  (see  earlier  references),  we  found  that  all  the  important  phenomena,  for  exam¬ 
ple  period  doubling  cascades  to  chaos,  are  completely  captured  by  a  truncated  system  -  this  was 
verified  by  comparison  of  numerical  solutions  of  the  infinite  dimensional  dynamical  system  (PDE)  and 
the  truncated  finite  dimensional  one,  with  emphasis  on  a  window  that  yields  a  sequence  of  eight  period 
doublings  before  chaos  sets  in.  From  hindsight  we  discovered  that  in  the  cases  of  the  smaller  values 
of  v  studied  here,  as  many  as  approximately  40  coupled  nonlinear  equations  must  be  retained  if  the 
time-periodic  solutions  are  to  be  described  with  satisfactory  accuracy. 

Of  particular  relevance  to  the  present  work  are  the  studies  of  Hyman  and  Nicolaenko  [18]  (also 
Hyman  et  al.  [19],  Kevrekides  et  al.  [22],  and  Jolly  et  al.  [20]).  These  studies  are  mostly  concerned 
with  the  extensive  numerical  study  of  the  related  KS  equation  for  the  integral  in  space  of  our  function 
u(x,t).  The  equation  used  is 

ui  +  4uxxxx  +  <*(««  +  \ux)  =  0  . 

(x,r)  e  R'xR+  ,  (1.2) 

u(x,0)  =  mo0O  ,  u{x+2n,t)  =  u(x,t)  . 

The  parameter  a  above  is  related  to  v  by  a  =  4/v.  The  range  of  a  used  was  4  £  a  £  120,  which 
corresponds  to  1/30  £  v  <,  1.  Various  attractors  were  constructed  by  changing  the  initial  condition  as 
well  as  a,  by  continuation  methods  where  the  initial  condition  is  renormalized  to  its  current  value 
before  a  is  changed.  The  bifurcation  diagram  reported  contains  windows  of  various  attractors  (fully 
modal,  bimodal,  trimodal  and  tetramodal  fixed  points)  separated  by  oscillatory  or  chaotic  orbits.  A 
time-periodic  orbit  is  reported  in  a  small  window  17.30  22.50. 

The  current  study  obtains  a  bifurcation  diagram  which  is  analogous  to  the  above,  but  some  essen¬ 
tial  differences  must  be  pointed  out.  First  of  all,  our  concern  is  not  with  the  full  description  of  the 
inertial  manifolds  which  are  by  no  means  easy  to  obtain  since  the  initial  conditions  belong  to  an 
infinite-dimensional  vector-space.  Instead  we  want  to  investigate  transition  to  chaos,  which,  as  we 
have  observed,  occurs  according  to  the  classical  routes  of  period-doubling  bifurcations.  Such 
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transitions  to  chaos  have  been  described  for  iteration  maps  (Feigenbaum  [10],  Collet  and  Eckman  [7], 
see  Devaney  [9]  for  a  full  description)  as  well  as  for  low  dimensional  systems  of  ordinary  differential 
equations  such  as  the  forced  Dufflng  equation  which  is  a  second  order  system,  and  for  the  Lorenz 
attractor  (see  Guckenheimer  and  Holmes  [16]  and  references  therein).  Numerical  evidence  of  a  period¬ 
doubling  sequence  has  been  presented  by  Knoblock  el  al.  [23]  and  Lennie  et  al.  [27]  for  the  two- 
dimensional  convection  in  a  Boussinesq  fluid.  Many  interesting  results  have  also  been  found  for  the 
Ginzburg-Landau  equation  by  Keefe  [21]  (see  also  Moon,  Huerre  and  Redekopp  [30],  Sirovich  and 
Rodriguez  [35])  including  a  period-doubling  sequence  that  supports  an  accurate  estimate  for  the 
Feigenbaum  number.  Finally  evidence  of  a  period-doubling  sequence  has  been  discovered  for  a  trun¬ 
cated  model  of  the  Kuramoto-Sivashinsky  equation  by  Brown  et  al.  [5]  (see  also  Jolly  et.  al  120]);  two 
period-doublings  are  reported  in  a  window  which  corresponds  to  the  tiny  range 
0.121223  ^  v  <.  0.121304  for  our  problem.  In  the  results  we  present  here  the  initial  condition  is  taken 
to  be  the  same  for  all  v.  Global  attractors  are  insensitive  to  such  a  choice,  but  in  regimes  where  the 
solution  becomes  time-periodic  and  period-doubling  bifurcations  take  place  the  present  choice  of  initial 
condition  enables  both  an  accurate  large-lime  integration  as  well  as  a  quantification  of  the  cascades 
(see  later  also).  We  emphasize  here  that  in  this  article  we  exclude  a  complete  study  of  the  phase-space 
of  initial  conditions. 

The  structure  of  the  paper  is  as  follows.  Section  2  describes  the  numerical  meUiods  used,  while 
in  Section  3  we  begin  to  describe  each  window  separately  by  showing  representative  numerical  experi¬ 
ments  in  each  case.  Here  we  describe  the  first  five  windows  which  support  solutions  that  are  steady  (i) 
(fully  modal,  bimodal,  fully  modal  again,  trimodal),  lime-periodic  (ii),  fully  modal  steady  (iii),  time- 
periodic  again  (iv)  and  steady  tetramodal  (v).  In  Section  4  results  for  the  sixth  window  and  beyond 
are  given;  the  sixth  window  is  time-periodic  and  contains  a  complete  infinite  sequence  of  period¬ 
doubling  bifurcations  leading  to  chaos.  Results  are  presented  for  both  a  direct  simulation  (classic?.! 
many  mode  discretization)  and  a  truncated  model  (traditional  Galerkin)  with  the  low  modal  beha  «or  of 
the  dynamical  system  clearly  brought  out.  The  results  are  used  to  compute  the  Feigenbaum  number  for 
the  cascade  which  is  found  to  be  in  excellent  agreement  with  Feigenbaum ’s  universal  theory.  This 
Section  also  includes  some  preliminary  results  of  solutions  at  smaller  v  showing  ,ie  appearance  of 
steady  attractors  of  decreasing  spatial  period;  lime-periodicity  and  chaotic  behavic.  persist  as  expected. 
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2.  Numerical  methods. 

As  mentioned  earlier,  our  concern  here  is  primarily  with  the  numerical  simulation  of  (1.1).  Two 
numerical  methodologies  are  used:  A  classical  numerical  discretization  with  many  degrees  of  freedom 
(which  we  term  the  "direct"  method)  and  a  traditional  Oalerkin  projection  onto  a  low-dimensional  iner¬ 
tial  manifold  that  yields  a  coupled  nonlinear  system  of  ODE’s  to  solve.  In  the  former  case  we  choose 
to  perform  the  time  integration  by  treating  the  linear  part  of  the  operator  in  Fourier  space,  and  the  non¬ 
linear  term  in  real  space  by  a  scheme  that  does  not  add  numerical  viscosity.  This  split  step  method 
maintains  second  order  accuracy  by  the  repeated  application  of  the  semi  groups,  as  suggested  by  the 
Strang  split  method.  The  time-step  is  chosen  adaptively  subject  to  several  controls  that  ensure  both 
numerical  stability  and  high  accuracy;  the  error  in  the  experiments  reported  here  is  kept  well  below  the 
10-6  level.  Such  requirement  is  necessary  to  capture  the  nonlinear  coupling  of  Fourier  modes  that  are 
undergoing  period-doubling  cascades.  For  the  Galcrkin  projection  the  Fourier  coefficients  were  com¬ 
puted  for  a  pre-determined  dimension;  a  fourth  order  Runge-Kutta  method  was  used  to  advance  the 
nonlinear  part  in  time,  while  the  linear  part  was  treated  exactly  by  use  of  the  split  method  described 
above  together  with  repeated  semi  group  application  that  maintains  overall  accuracy.  A  natural  accu¬ 
racy  test,  throughout  the  computation,  is  the  conservation  of  the  spatial  integral  of  u .  The  initial  condi¬ 
tion  used  is 

u0(x)  =  -sinU)  , 

and  a  consequence  of  this  is  that  the  integral  over  a  spatial  period  vanishes  for  all  time,  and  further 
more  the  solution  maintains  its  odd  parity.  More  generally,  if  the  spatial  profile  of  the  solution  is  ini¬ 
tially  of  the  form 

u  (*  ,0)  =  £  at  sin  kx  , 

*21 

then  the  mechanisms  of  the  evolution  of  the  partial  differential  equation  respect  this  structure  since, 

ua  +  v  {k2  -  vk4)  ak  sin  kx 

*21 

and  the  nonlinear  part  becomes, 

=  (y«2),  =  y-y-(£  <**  sin  kxf  = 

*21 

=  •  •  •  =  -  4 -  7  E  sin  kx. 

*21  X21  x+n^* 

Thus  the  spatial  profile  remains  a  linear  combination  of  sines  for  all  times,  i.e. 
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u(x,t)  =  £  ak(t )  sin  kx  . 

kil 

Due  to  the  above  observation  the  PDE  (1.1)  becomes  equivalent  to  an  infinite  dimensional  system  of 
ordinary  differential  equations  for  ak's,  i.e. 

ak  =  ( k 2  -  v*4)  ak  +  ~  \  X  for  k  £  1  .  (2.1) 

tel  ten-* 

The  odd-parity  property  was  utilized  by  a  suitable  filtering  of  Fourier  coefficients;  the  filtering  toler¬ 
ance  was  set  to  a  value  between  KT10  and  10-12.  In  the  case  of  the  direct  computation  the  spatial 
discretizations  were  done  on  grids  ranging  from  128  to  1024  points.  For  windows  other  than  the  first 
steady  attractor,  at  least  512  points  were  used.  Results  that  are  obtained  by  the  projection  method  will 
be  pointed  out  together  with  the  number  of  modes  retained  for  the  dynamical  system. 

The  results  presented  here  were  produced  after  several  hundreds  of  numerical  experiments  which 
tracked  the  time  evolution  of  the  profile.  Us  Fourier  modes,  its  energy  (L2-norm)  and  the  phase-plane 
of  the  energy  and  in  some  marginal  cases  of  certain  Fourier  coefficients.  The  lime  evolution  of  the 
energy  is  very  representative  of  the  structure  of  the  various  attractors.  Steady-state  attractors  evolve  to 
a  constant  large-time  state,  while  periodic  or  chaotic  solutions  have  a  time-oscillatory  energy  profile. 
Near  window  and  subwindow  (period-doubling)  boundaries  for  the  smaller  values  of  v  at  least,  the 
time  integration  required  as  many  as  200  time  units  (8xl05  time  steps)  in  order  to  sharply  classify  the 
v-windows  and  subwindows.  All  periodic  solutions  reported  here  were  derived  by  a  phase-plane 
analysis  of  the  data.  As  is  described  in  the  next  Section  many  period-doubling  bifurcations  are  evi¬ 
denced  from  the  phase-plane,  which  is  used  in  conjunction  with  energy  plots  to  build  the  geometrical 
structure  of  the  solution  in  one  time-period  and  so  accurately  estimate  the  the  v-posilion  of  the  hyper¬ 
bolic  point  corresponding  to  a  period-doubling  bifurcation,  for  example.  Table  1  presents  a  summary 
of  all  windows  discovered  so  far,  together  with  sharp  computed  estimates  of  window  boundaries. 

3.  Detailed  numerical  experiments.  Pre-chaotic  windows. 

When  the  parameter  v  decreases  below  unity,  the  solution  undergoes  its  first  Hopf  bifurcation.  In 
this  Section  we  give  a  description  of  the  dynamics  in  each  of  the  first  five  windows. 

(i)  Steady  attractors  ,  .05985  <,  v  <  1. 

The  main  characteristic  of  this  window,  is  the  attraction  of  the  initial  profile  to  a  steady-state.  The 
steady-state  is  approached  exponentially  fast  and  is  numerically  achieved  after  a  relatively  short  time 
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(5-10  time  units  in  most  cases).  However  the  numerical  convergent  to  the  steady  steady  state  is 
significantly  delayed  towards  the  lower  end  of  the  window  where  the  first  periodic  attractor  sets  in.  We 
subdivide  the  window  into  the  following  four  subwindows  : 

(a)  A  fully  modal  steady  attractor  when  .248  <.  v  <  1  . 

(b)  A  bimodal  steady  attractor  (l.e.  where  all  the  odd  frequencies  vanish)  when  .0758  ivS  .247  . 

(c)  Another  fully  modal  steady  attractor  when  .06697  Svi  .0758  . 

(d)  A  trimodal  steady  attractor  (where  only  the  frequencies  which  are  multiples  of  three  survive)  when 
.05985  <;  v  £  .06695  . 

For  values  of  v  near  the  lower  boundary  in  fact,  we  find  that  there  are  about  30  modes  which 
have  an  amplitude  larger  than  KT6.  In  Figures  (1.1),  (1.2)  and  (1.3)  we  give  sets  of  these  attracting 
profiles  for  different  decreasing  values  of  v;  Figures  (1.1a)  and  (1.1b)  cover  ranges  of  v  in  the  first 
subwindow.  Figure  (1.2a)  covers  a  range  of  v  at  the  lower  end  of  the  first  subwindow  as  well  as  at 
the  upper  end  of  the  second  subwindow  indicating  a  smooth  transition  from  the  fully  modal  first 
subwindow  to  the  bimodal  second  subwindow.  Figure  (1.2b)  also  covers  a  range  of  v  in  the  second 
subwindow,  whereas  figure  (1.3a)  covers  a  range  in  the  third  subwindow  where  the  attracting  profiles 
are  trimodal.  Comparison  of  figures  (1.2b)  and  (1.3a)  provides  the  indication  that  the  transition  of  the 
attracting  profiles  between  tire  two  subwindo.vs  is  smooth.  Finally  figure  (1.3b)  covers  a  range  of 
values  of  v  in  the  fourth  and  last  subwindow  where  the  attracting  steady  profiles  are  trimodal.  The 
change  in  the  shape  of  the  limiting  profiles  in  this  subwindow  is  hardly  visible,  even  though  v  changes 
"significantly"  from  .06695  to  .06  .  This  is  not  surprising  since  under  a  similarity  rescaling  (see  below) 
they  correspond  to  fully  modal  steady  profiles  for  v  between  .54  and  .603  where  not  much  change  Is 
observed  either. 

For  values  of  v  which  are  near  1,  the  profile  is  qualitatively  similar  to  a  sine  wave.  As  v 
decreases,  however,  one  more  maximum  and  one  more  minimum  are  obtained  after  the  profile  attains 
an  inflection  at  the  anti-symmetry  line  x  =  n,  Such  solutions  were  also  obtained  by  Frisch  et  a 1.  [14] 
for  values  of  v  larger  than  or  equal  to  0.35;  as  a  result  of  this  they  did  not  see  profiles  with  four  turn¬ 
ing  points.  It  is  interesting  at  this  point  to  describe  an  important  property  of  the  steady  KS  equation;  If 
u(x )  is  a  solution  corresponding  to  a  viscosity  v  of  the  ODE  obtained  by  dropping  the  time  depen¬ 
dence  in  (1.1),  then  for  any  integer  N  there  exists  another  solution  Nu{Nx)  with  viscosity  v/A2.  Due  to 
the  periodic  boundary  conditions  used  here,  N  must  be  an  integer.  Tills  self-similar  property  of  the 
ordinary  differential  equation,  satisfied  by  the  steady  states  of  the  KS  equation,  permits  the  generation 
of  a  one  parameter  family  of  steady-states  which  grow  in  amplitude  and  decrease  in  period.  In  sub- 


section  (v),  the  tetramodal  attractor,  we  give  a  numerical  example  of  this  self-similar  property  connect¬ 
ing  steady-stales  for  N  =  4. 

(ii)  Time  periodic  attractors.  .055238  £  v  <,  .0599. 

The  large  time  evolution  of  our  initial  value  problem  for  v  lying  in  the  range  covered  by  this 
window,  is  attraction  to  a  time-periodic  orbit.  In  all  the  results  presented  here,  the  period  is  defined  for 
the  energy  (  L2-norm  )  of  the  solution.  It  is  clear  that  the  solution  itself  can  have  half  the  period  of 
the  energy  but  such  a  case  was  not  observed  when  the  Fourier  coefficients  were  considered  individu¬ 
ally.  Periodicity  is  determined  by  examination  of  the  energy  phase-plane  (  E  versus  dEldt  ),  and  the 
period  is  evaluated  from  the  evolution  of  the  energy  with  time.  All  the  solutions  in  this  window  are 
fully  modal,  and  are  characterized  by  a  time-periodicity  (the  time  period  will  be  denoted  by  x  )  which 
begins  at  a  value  of  approximately  x  =  1.0183,  Increases  slowly  to  x  =  1.18  at  a  value  of  v  =  .059  and 
then  doubles  to  x  =  2.374  at  a  value  of  v  =  .0559.  The  doubling  Is  instantaneous,  in  parameter  space, 
and  comes  about  by  a  slight  disalignment  of  the  maxima  and  minima  of  the  energy-time  graph.  Some 
representative  results  from  this  window  are  given  in  Figures  (2.1a,b)-(2.3a,b).  The  solution  undergoes  a 
period  doubling  as  is  evidenced  by  comparison  of  the  phase-planes  at  v  =  057  and  v  =  .0559  -  see 
Figures  (2.1a),  (2.2a).  Figures  (2.3a,b)  correspond  to  v  =  .0555  and  it  can  be  seen  that  the  4  cycle  pro¬ 
perty  of  the  phase  diagram  is  much  more  pronounced  than  at  the  slightly  higher  value  of  v  =  .0559; 
such  drastic  changes  in  the  dynamics  over  small  domains  of  the  parameter  space  are  typical  and  are 
encountered  later  on  also.  A  summary  of  the  results  from  this  window  is  contained  in  Table  2. 

(iii)  Steady  fully  modal  attractors.  .03965  <,v  Z  .055235. 

Here  a  new  steady  window  is  found  with  fully  modal  solutions  which  do  not  correspond  to  a  res¬ 
caled  value  of  the  parameter  v  under  the  self-similar  scaling  described  earlier.  The  new  element  here 
is  the  extremely  slow  convergence,  still  exponential  nevertheless,  towards  the  fixed  point  attractor.  At 
such  values  of  v  the  equations  must  be  integrated  to  long  times;  this  is  to  be  expected  since  a  straight¬ 
forward  perturbation  analysis  above  the  window-boundary  yields  a  rale  of  decay  of  the  distance  (for 
example  die  L2-norm)  from  the  attractor,  which  is  inversely  proportional  to  the  perturbation  of  v  from 
the  boundary.  This  is  illustrated  in  Figures  (3.1a)-(3.1c)  for  v  =  .0397.  Figure  (3.1a)  shows  the  varia¬ 
tion  of  die  energy  with  time,  while  Figures  (3.1b)  and  (3.1c)  show  die  variation  with  time  of  the  L°° 
norm  and  the  L~  norm  of  the  gradient  of  the  solution  respectively.  It  can  be  seen,  that  Initially  the 
soludon  passes  through  an  oscillatory  and  almost  periodic  stage,  and  at  a  time  of  approximately  20 
units  its  amplitude  follows  an  envelope  that  takes  the  profile  into  the  stationary  attractor.  In  fact 
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during  the  oscillatory  stage,  the  phase-plane  of  the  energy  is  similar  to  that  of  the  periodic  window 
described  next;  once  the  steady-state  begins  its  attraction  the  phase  curve  converges  to  a  sink  of  course. 

(iv)  Time  periodic  attractors.  ,03735  £  v  £  .0396. 

Periodicity  is  again  determined  by  examination  of  the  energy  phase-plane,  and  likewise  the  period 
is  computed  from  the  energy  versus  time  plots.  In  this  window  the  period  begins  at  a  value  of  approx¬ 
imately  x  =  .515  when  v  =  .0396,  with  a  very  slow  convergence  to  the  periodic  attractor,  and  it 
increases  smoothly,  in  parameter  space,  until  it  reaches  the  value  x  =  .684  when  v  =  .03798  with  a 
much  faster  convergence  to  the  the  limit  circle;  the  period  then  suddenly  doubles  to  t  =  1.37  at  a  value 
of  v  =  .03795.  As  v  is  decreased  further,  the  period  undergoes  a  new  slow  increase  to  reach  a  value 
x  =  1.43  at  v  =  .0377;  a  small  decrease  of  v  to  a  value  .0376  gives  a  new  period  doubling  with 
x  =  2.87.  The  period  again  increases  slowly  to  a  value  x  =  3.0  at  v  =  .03736.  At  this  point  a  slight 
decrease  of  v  to  the  value  v  =  .037355  gives  a  period-halving  with  the  energy  having  a  period  x  -  1.5. 
The  last  value  of  v  computed  in  this  window  is  v  =  .03735,  for  which  the  period  is  x  =  1.5.  The  next 
numerical  experiment  with  v  =  .037345  is  outside  this  time-periodic  window  and  yields  a  tetramodal 
steady-state  (see  later).  These  results  are  summarized  in  Table  3. 

The  results  of  Table  3  were  obtained  after  numerous  computations  together  with  a  careful  exami¬ 
nation  of  the  oscillatory  dynamics  by  means  of  phase-plane  analysis,  and  representative  solutions  are 
given  here.  Some  cases  are  described  in  full  including  the  variation  and  phase-planes  of  certain  Fourier 
modes,  while  other  cases  are  described  by  the  energy  time  series  amid  the  energy  phase-plane.  These 
two  plots  are  enough  to  conclude  periodicity  of  solutions  and  their  periods.  In  Figures  (4.1a-d)  we 
present  results  for  v  =  .038  which  is  a  representative  parameter  at  the  beginning  of  the  window.  These 
Figures  show  the  variation  of  the  energy  with  time,  the  energy  phase-plane,  and  the  phase  planes  of  the 
first  and  tenth  Fourier  modes  respectively;  time  periodicity  or  solutions  is  clear  and  we  note  that  the 
Energy  phase  plane  (Figure  (4.1b))  is  traced  out  21  times  since  the  time  interval  it  depicts  contains  21 
periods.  The  time  period  of  Fourier  coefficients  is  the  same  as  that  for  the  energy,  a  fact  that  need  not 
hold  in  general  (see  earlier  comments);  the  reason  it  holds  in  the  cases  presented  by  this  article,  is  that 
the  actual  solution  u(x,t)  has  a  non-zero  mean  part  in  general. 

As  has  been  mentioned  already  (Table  3),  a  period-doubling  takes  place  at  v  =  .03795.  This  is 
evidenced  very  clearly  in  the  phase-plane  of  the  energy  for  this  case  which  is  given  by  Figure  (4.2a). 
The  phase-plane  consists  of  two  loops  which  are  initially  almost  coincident  (in  fact  they  were  coin¬ 
cident  at  v  =  .038)  and  this  change  produces  the  period-doubling.  The  corresponding  energy  variation 
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wiih  time  can  be  found  in  Figure  (4.2b),  from  which  the  period  is  calculated  to  have  a  value  of  1.37. 
The  mechanism  of  the  period  doubling  can  also  be  seen  from  (4.2b),  where  every  other  wave  crest  and 
trough  are  reduced/increased  slightly  thus  producing  a  phase-plane  of  index  two  and  a  period-doubling. 
We  note  also  that  there  is  no  abrupt  change  in  the  wave-form  that  undergoes  a  period-doubling;  only 
the  value  of  the  period  doubles  abruptly.  Next  we  reduce  the  value  of  v  by  5x1  O'5  to  v  =  .0379.  We 
are  still  in  the  same  subwindow,  but  the  phase  plane  diagram  is  now  much  more  pronounced  as  Figure 
(4.3a)  shows.  The  results  for  the  last  v  from  this  subwindow  is  presented  in  Figures  (4.3b),  (4.3c)  and 
(4.3d)  which  correspond  to  the  phase-planes  of  the  energy,  the  first  and  sixth  Fourier  coefficients 
respectively,  at  a  value  of  v  =  .0377. 

A  reduction  of  v  to  the  value  .0376  brings  us  into  a  new  subwindow  with  the  solution  undergo¬ 
ing  a  period  doubling.  Evidence  for  this  can  be  found  in  Figure  (4.4a)  which  depicts  the  phase-plane 
of  the  Energy  at  v  =  .0376;  the  time  period  is  t  =  2.87;  Figures  (4.4b)  and  (4.4c)  are  also  given  to 
show  the  evolution  of  the  sixth  Fourier  coefficient  with  time  and  its  corresponding  phase-plane.  The 
lime-periodicity  is  again  clear  and  similar  results  were  found  for  all  leading  order  (i.e.  large  enough  for 
graphical  purposes)  Fourier  coefficients.  Such  results  are  of  course  expected  but  they  serve  as  a  strict 
check  on  the  numerical  accuracy  of  the  solutions  so  that  we  can  confidently  state  that  the  time  oscilla¬ 
tions  are  by  no  means  spurious,  but  in  fact  persist  in  a  regular  and  accurate  manner.  The  last  case 
computed  in  this  subwindow  has  v  =  .03736.  The  energy  evolution  and  phase  plane  curves  are  given 
in  Figures  (4.5a)  and  (4.5b).  The  solution  is  just  about  to  undergo  a  period  halving  and  the  index  two 
property  of  the  phase  plane  is  only  just  discernible.  A  reduction  of  v  to  the  value  .037355  gave  a 
period  halving,  and  Figure  (4.5c)  shows  the  energy  phase-plane  for  the  last  member  of  this  window 
that  was  computed,  v  =  .03735.  .  As  expected,  the  mechanism  that  produced  period-doublings  is  now 
working  the  opposite  way  to  produce  period-halvings  (similar  windows  were  also  found  by  Knoblock 
et  al.  [23]  for  a  convection  problem).  As  mentioned  earlier,  a  numerical  experiment  with  v  =  .037345 
produced  a  tetramodal  steady-state  which  is  representative  of  the  window  described  next. 

It  is  also  interesting  to  observe  the  spatio-temporal  evolution  of  the  profiles,  and  we  present  this 
for  two  representative  values  of  v  from  this  window,  v  =  .039  and  .0375  respectively.  The  way  this  Is 
done  is  by  plotting  on  the  same  diagram  the  solution  at  successive  increasing  limes.  To  make  the  evo¬ 
lution  clearer  the  profiles  are  shifted  vertically  at  each  time  interval  producing  pictures  that  capture  the 
spatio-temporal  evolution  and  in  particular  the  time-periodicity  of  the  solutions.  Figure  (4.6) 
corresponds  to  v  =  .039  while  Figure  (4.7)  shows  the  results  for  v  =  .0375. 


-10- 


(v)  Steady  tetramodal  attrac' -  -s.  .037348  £  v  $  .0344. 

We  present  only  one  case  (v  =  .0373)  from  this  window  since  the  results  are  qualitatively  similar, 
the  difference  being  in  the  final  amplitude  of  the  solution.  For  v  =  .0373,  therefore,  a  steady-state  Is 
reached  after  approximately  seven  time  units.  The  evolution  of  the  energy  is  shown  in  Figure  (5.1a) 
from  which  it  is  seen  that  after  a  few  oscillations  in  time,  the  solution  gets  attracted  to  a  tetramodal 
steady-state;  this  was  determined  from  inspection  of  the  Fourier  coefficients.  The  corresponding  profile 
at  t  -  10  is  depicted  in  Figure  (5.1b)  which  clearly  shows  die  tetramodal  character  of  the  solution. 
This  means  that  die  solution  is  spatially  periodic  and  of  period  n/2.  Further  more,  this  solution  is 
obtainable  from  the  self-similar  re-scaling  (see  earlier)  of  the  bimodal  steady  stale  at  v  =  .’  '92  with 
die  value  of  N  -2.  This  was  checked  and  agreement  is  excellent. 

4.  Complete  period-doubling  cascade  to  chaos  and  computation  of  the  Feigenbaum  number. 

(vi)  Time  periodic  attractors.  Complete  period-doubling  sequence.  .0299756  <*\  <,  .0343. 

This  window  is  to  our  view,  the  most  important  found  in  this  work,  and  for  this  reason  we  con¬ 
centrate  on  die  results  it  produces  rather  than  try  to  sharpen  and  classify  the  bifurcation  diagram  at 
higher  values  of  the  viscosity  (the  classification  of  die  bifurcation  diagram  at  high  values  of  v  has  been 
the  subject  of  extensive  numerical  and  analytical  studies  by  various  authors  -  [18], [19],  [20],[2 1  ]).  The 
importance  of  diis  window  is  two  fold.  First,  it  supports  solutions  that  start  as  laminar  steady -states 
(tetramodal  waves)  and  dicn  become  time-periodic  (this  is  also  die  case  for  window  (iv)  above). 
Second,  the  time-periods  of  the  soludon  undergo  a  complete  sequence  of  period-doubling  bifurcations 
as  the  viscosity  v  is  decreased,  resulting  in  an  accumulation  point  beyond  which  the  solution  is  chaodc. 
Remarkably  enough  as  many  as  5  period -doubling  bifurcations  were  accurately  computed  by  the  direct 
method  and  8  period-doublings  for  the  Galerkin  procedure.  This  enables,  for  the  first  time,  the  accu¬ 
rate  computation  of  the  Feigenbaum  number  for  the  Kuramoto-Sivashinsky  equation.  Other  investiga¬ 
tions  of  PDE  systems  found  results  that  were  consistent  with  the  universal  Feigenbaum  constant  but  no 
conclusive  evidence  was  shown  (see  for  example  Knoblock  et  al.  [22]  for  a  convecdon  problem,  and 
Jolly  et  al.  [20]  for  the  KS  equation).  Significant  results  were  obtained  by  Keefe  [21],  however,  for 
the  Ginzburg-Landau  equation.  After  extensive  numerical  computations  of  even  periodic  solutions 
Keefe  found  a  window  containing  a  complete  sequence  of  period-doubling  bifurcations;  the  first  five 
bifurcations  were  accurately  obtained  and  based  on  the  results  the  Feigenbaum  number  for  the  cascade 
was  estimated  to  be  approximately  4.4  to  4.7.  The  numerical  method  u.  d  in  [21]  was  a  pseudo- 
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spectral  one  with  a  32  mode  discretization.  We  mention  here  that  these  bifurcations  take  place  on  very 
small  scales,  in  parameter  space,  and  require  a  lot  of  care  numerically,  as  well  as  an  extensive  compu¬ 
tational  effort  towards  a  fairly  sharp  estimation  of  subwindow  lengths  (a  subwindow  should  be  under¬ 
stood  to  be  a  region  In  parameter  space  that  separates  period-doubling  parameter  values).  We  begin  tills 
section  by  presenting  the  results  of  the  direct  method  by  use  of  at  least  512  Fourier  modes.  Most  of  the 
evidence  from  this  approach  is  graphical  and  serves  to  show  numerous  period  doublings  through 
phase-plane  analysis.  We  also  re  computed  this  same  window  by  means  of  a  truncated  model  that 
retains  the  leading  36  Fourier  modes.  The  results  of  the  two  methods  agree  to  at  least  4  decimal 
places  (by  results  here  we  mean  the  positions  in  v-space  of  the  various  subwindows);  a  quantitative 
analysis  of  the  results  obtained  from  the  36-mode  Galerkin  expansion  enables  the  accurate  computation 
of  tiie  universal  Feigenbaum  number  of  the  cascade  (see  later  for  details). 

The  direct  computation  of  the  PDE  provides  results  for  the  first  six  subwindows  over  which  the 
solution  attains  five  period  doublings.  Improvement  to  a  sixth  period  doubling  would  require  a  substan¬ 
tial  amount  of  computational  resources  and  is  wise  at  this  point  to  turn  to  finite  dimensional  versions 
of  the  equations.  For  example  computations  that  retain  36  Fourier  modes  are  in  general  at  least  8  times 
faster  than  direct  computations;  in  the  absence  of  sharp. analytical  remits  direct  computations  are  essen¬ 
tial  in  gauging  the  accuracy  or  sufficiency  of  truncated  models.  As  mentioned  already,  36  modes  are 
capable  of  capturing  a  complete  set  of  period-doubling  bifurcations  that  precede  chaos,  and  more 
importantly  are  in  complete  agreement  with  the  Feigenbaum  universal  constant.  To  give  an  example  of 
the  numerical  restrictions  that  are  present  in  either  direct  or  finite  dimensional  approaches,  our  direct 
numerical  solutions  for  v  =  .0299756  gave  a  time  periodic  solution  with  x  -  7.04  while  at  v  =  .02997 
the  solutions  appear  to  be  chaotic.  In  a  distance  of  less  than  5.6x10  5,  therefore,  numerous  period¬ 
doubling  bifurcations  take  place  (to  conjecture  an  infinite  number  would  be  too  hasty  because  the 
chaotic  solution  at  v  =  .02997  can  be  thought  of  as  an  oscillatory  solution  with  a  very  large  period;  the 
period  Is  so  large  that  discourages  any  attempt  at  its  numerical  computation).  It  can  be  seen,  therefore, 
that  we  would  very  soon  run  out  of  computational  accuracy,  since  subwindow  lengths  would  be  of  the 
same  order  as  computer  round-off  error.  We  decided,  therefore,  to  concentrate  on  the  first  five  subwin¬ 
dows  (eight  for  the  36  mode  model)  and  compute  the  Feigenbaum  number  from  these. 

For  brevity  of  the  description  we  will  restrict  our  presentation  to  phase-plane  curves  of  the 
energy.  Each  case  presented  was  computed  with  512  Fourier  modes  and  approximately  4x  1  ()4  time 
steps  (a  run  takes  about  20  minutes  CPU  time  on  a  CRAY-11  and  about  4  hours  CPU  time  on  a 
SPARC  work  station;  all  runs  use  double  precision).  Accuracy  tests  were  also  made  for  many  runs  by 
use  of  1024  Fourier  modes;  all  results  reported  here  are  insensitive  to  such  refinements  as  tire  truncated 
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system  also  shows  (see  later).  We  give  representative  results  from  each  subwindow  and  show  the 
behavior  near  subwindow  boundaries.  The  steady  tetramodal  solution  gets  attracted  to  a  time-periodic 
orbit  when  v  is  reduced  below  a  value  of  .0345.  Hie  phase-plane  of  the  L2-norm  of  the  solution  is  a 
1 -cycle  initially.  In  Figure  (6.1a)  the  phase  plane  for  v  =  .0342  is  shown  over  a  time  interval  from  1.0 
to  34.5.  It  is  clear  that  the  solution  gets  attracted  to  a  time-periodic  orbit  of  index  1;  the  period  is  0.37. 
As  v  is  decreased  the  period  increases  smoothly  and  the  1 -cycle  topology  of  the  phase-plane  persists 
until  the  first  period  doubling  suddenly  appears  at  a  value  of  v  just  above  .0304.  The  phase-plane 
corresponding  to  v  =  .0304  is  not  shown  because  the  period-doubling  bifurcation  has  just  taken  place 
and  within  graphical  accuracy  the  phase-plane  appears  to  have  the  same  topology  as  in  the  previous 
subwindow.  Careful  inspection  of  the  left  side  of  the  curve,  however,  shows  that  there  are  two  distinct 
branches;  the  period  of  the  solution  is  0.89.  The  two  branches  begin  to  move  apart  as  v  decreases 
further.  An  example  of  this  change  in  the  phase-plane,  which  in  turn  corresponds  to  a  change  in  the 
oscillatory  structure  of  the  solution,  is  seen  in  Figures  (6.1b-c);  the  value  of  v  is  .0303  and  the  period 
is  the  same  as  before,  x  =  0.89.  The  last  picture  we  present  from  this  subwindow  is  Figure  (6. Id)  for  a 
value  v  =  .03007.  The  characteristics  of  the  time  oscillations  are  qualitatively  the  same  as  before  with 
the  two  loops  well  separated  and  with  time-period  x  =  0.88. 

We  next  decrease  the  value  of  v  slightly.  A  new  subwindow  is  found  which  develops  after  a 
period-doubling  bifurcation  at  a  value  of  v  between  .03007  and  .030065.  The  phase-plane  curve  now 
becomes  a  4-cycle  and  in  one  complete  time  oscillation  the  energy  has  eight  turning  points,  four  of 
which  are  maxima  and  four  minima.  A  characteristic  phase-plane  curve  for  bifurcation  parameters  near 
the  upper  edge  of  the  second  subwindow,  is  Figure  (6.2a)  which  has  v  =  .03006;  the  time-period  is 
x  =  1.76.  The  index  4  property  is  evident  and  appears  by  a  slight  disalignment  of  the  phase  plane 
branches.  We  also  note  that  at  this  value  of  v  the  maxima  of  the  energy  are  noticeably  different 
whereas  the  minima  are  almost  indistinguishable.  Other  cases  from  this  subwindow  are  also  given  in 
Figures  (6.2b-d).  The  value  of  the  viscosity  is  v  =  .030015  ,  .03  ,  .029998  respectively.  The  period  is 
1.76  and  the  4-cycle  property,  characteristic  of  this  subwindow,  is  very  clear.  Figure  (6. 2d),  which  has 
v  =  .029998  is  at  the  lower  edge  of  the  second  subwindow  and  Is  representative  of  the  dynamics  just 
before  a  second  period-doubling  bifurcation. 

Next  the  value  of  v  was  decreased  to  .0299975;  at  first  sight  this  v  seems  to  belong  to  the  second 
subwindow  but  a  careful  inspection  of  one  of  the  loops  shows  that  a  period-doubling  takes  place  some¬ 
where  between  v  =  .029998  and  .0299975  leading  to  an  estimated  time-period  of  x  =  3.52  and  phase- 
plane  characteristics  for  this  subwindow  of  index  8.  For  smaller  values  of  v  the  phase-plane  loops 
begin  to  drift  apart.  Figures  (6.3a-c)  show  a  sequence  of  energy  phase-plane  diagrams  for  three 
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decreasing  values  of  v  that  belong  to  this  third  subwindow;  these  values  are 
v  =  .029995  ,  .02999  ,  .0299815.  The  last  computation  in  this  subwindow  is  at  v  =  .0299815  and  the 
period  is  x  =  3.52. 

The  next  numerical  experiment  has  v  =  .029981.  Once  more,  careful  inspection  of  the  phase- 
plane  shows  that  a  16-cycle  is  obtained  and  hence  another  period-doubling  bifurcation  takes  place 
somewhere  between  v  =  .0299815  and  .029981.  To  sc<  this  one  must  concentrate,  for  example,  on  that 
region  of  the  phase  plane  described  by  17.5  <  E  <  18  and  dEldt  near  zero.  Before  the  period-doubling 
bifurcation  takes  place  there  is  a  single  branch  traced  out  in  the  phase-plane  while  now  two  branches 
are  distinguishable  within  the  limits  of  graphical  accuracy.  At  smaller  values  of  v  the  branches  move 
apart  and  the  16-cycle  property  is  very  clearly  seen.  Figures  (6.4a-c)  give  representative  phase-planes 
corresponding  to  a  sequence  of  viscosities  which  span  this  fifth  subwindow.  The  bifurcation  parameters 
are  v  =  .02998  ,  .029978  ,  .0299756  respectively.  The  value  v  =  .0299756  is  also  the  last  computation 
from  this  subwindow;  the  period  is  x  =  7.04. 

In  Figures  (6.5a-d)  we  present  results  after  another  period-doubling  occurs.  The  first  case,  Figure 
(6.5a),  has  v  =  .029975  and  the  phase-plane  appears  to  be  a  32-cycle  with  a  period  of  at  least  14.0. 
The  complexity  of  the  dynamics  becomes  more  enhanced  as  v  is  decreased  even  further,  and  it 
becomes  very  difficult  and  computationally  expetisive  to  try  and  follow  more  period  doubling  bifurca¬ 
tions  than  the  ones  described  here.  Figures  (6.5b-d)  give  a  flavor  of  the  complicated  dynamics  that  one 
has  to  resolve,  for  the  values  v  =  .02997  ,  .02995  ,  .0299  respectively. 

Computation  of  the  Felnenbaum  number. 

As  mentioned  earlier,  the  results  we  have  front  the  direct  numerical  simulation  are  probably  not 
extensive  enough  for  an  accurate  determination  of  the  Feigenbaum  number.  The  lacking  ingredient  is 
the  sharp  estimation  of  subwindow  boundaries  which  take  on  a  central  role  in  such  a  compulation.  As 
a  result  of  this,  a  straightforward  calculation  of  the  ratio  of  successive  subwindow  lengths  suffers  from 
certain  inaccuracies.  For  completeness,  however,  c  give  the  results  and  stress  that  they  can  be 
Improved  by  resort  to  a  much  more  extensive  numerical  search.  If  we  label  the  lengths  of  successive 
subwindows  by  /;  ,  /= 1,°°  then  the  following  ratios  are  obtained  from  the  results  of  Table  4  ; 
/2// 1  =  5.42,  l3ll2  =  1135  and  ljl3  =  2.77.  The  cascade  is  still  in  its  early  stages  but  the  distance,  in 
v-space,  of  the  fourth  subwindow  lower  boundary  fiom  the  point  where  the  strange  attractor  seems  to 
commence,  is  of  the  order  of  1IT6.  Due  to  the  non-sharp  classification  of  the  subwindows,  the  ratios 
are  only  of  the  same  order  of  magnitude  as  the  universal  Feigenbaum  constant.  Our  data  is  much  more 
extensive  for  the  36-mode  Galerkin  approach  and  they  yield  excellent  results  as  we  explain  next. 
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The  period-doubling  route  to  chaos  is  well  documented  for  certain  low  dimensional  models  (see 
earlier  references),  but  it  is  believed  that  the  Feigenbaum  constant,  originally  found  for  iteration  maps, 
is  universal  and  appears  in  other  high-dimensional  systems  also.  In  what  follows  we  present  numerical 
evidence  of  tills  universality  for  the  Kuramoto-Sivashinsky  equation.  Suppose  that  period-doublings 
take  place  at  some  successive  values,  v„  say,  of  the  viscosity  parameter.  The  theory  predicts  an  accu¬ 
mulation  of  such  bifurcation  points  as  n  — >  °°  and  in  fact  the  universal  constant  is  given  by  : 

lim  V--~- V—  =  4.6692016...  . 

"  Vn  +  1  _  Vn 

V.  > 

The  obvious  way  to  test  our  results  against  the  theory  is  to  form  the  successive  ratios  indicated  above, 
for  our  computed  subwindow  boundaries.  The  results  of  this  approach  are  given  in  Table  5.  The  first 
column  of  the  Table  gives  our  estimate  for  subwindow  boundaries  (for  example  the  first  two  numbers 
in  the  first  column  show  the  beginning  and  end  respectively  of  the  first  periodic  subwindow;  the 
second  and  third  numbers  indicate  the  beginning  and  end  of  the  second  periodic  subwindow  etc.),  the 
second  column  shows  the  subwindow  length  and  the  third  column  contains  the  ratio  of  successive 
subwindow  lengths.  Note  that  the  cascade  develops  as  we  move  down  the  Table  and  by  Uie  time  the 
third  period-doubling  bifurcation  is  attained  the  geometric  ratio  of  successive  bifurcation  points  is  4.27, 
a  number  already  close  to  Feigenbaum’s  universal  constant.  It  is  worth  noting  that  subwindow  lengths 
are  shrinking  significantly  and  the  last  subwindow  that  was  computed  accurately  had  a  length  of 
1.3x10  .  It  is  clear,  therefore,  that  the  difficulty  of  the  numerical  experiments  increases  as  the  cascade 
develops  and  some  special  continuation  algorithms  were  used  towards  the  expediency  of  our  estimates. 
Such  accuracy  measures  are  essential  in  the  production  of  four  successive  bifurcation  ratios  from  the 
fourth  to  tire  eighth  subwindow,  which  are  at  an  error  of  at  most  4%  from  tire  universal  constant  of 
4.669201609...  .  Various  other  tests  can  also  be  performed  using  the  data.  For  example  we  can  pick  a 
subwindow  which  begins  at  v„  and  which  has  length  /„ .  Then  by  use  of  the  computed  accumulation 
point  vc  =  .029969275  and  by  assumption  of  an  infinite  cascade,  it  is  easy  to  show  that  there  exists  a 
geometric  ratio 


which  takes  the  point  v„  to  the  accumulation  point  vc.  The  number  8  above,  provides  then  the 
"Feigenbaum  number"  for  the  chosen  subwindow.  The  theory  tells  us  that  as  the  value  of  v„  is  chosen 
to  be  progressively  closer  to  the  accumulation  point,  8  then  converges  to  the  universal  constant.  This 
version  of  the  cascade  theory  was  also  applied  to  our  numerical  data.  The  results  are  similar  to  those  in 
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Table  5  but  they  are  slightly  more  accurate,  the  gain  in  accuracy  is  due  to  elimination  of  combination 
of  errors  that  a  successive  ratio  inevitably  contains,  Finally,  we  also  computed  the  Feigenbaum  number, 
based  on  subwindows  which  are  not  successive  but  have  a  factor  of  4  difference  in  their  time-periods. 
A  ratio  is  again  obtained  whose  square  root  should  converge  to  the  universal  constant  after  an  infinite 
number  of  bifurcations.  After  eight  bifurcation  (even  less  in  fact)  the  numbers  obtained  are  again  very 
similar  to  those  of  Table  5  and  we  omit  further  details  here. 

Finally  in  this  Section  we  present  a  brief  summary  of  some  of  the  results  obtained  for  values  of  v 
which  take  the  dynamics  beyond  the  compiete  sequence  of  period-doubling  bifurcations  just  described. 
The  richness  of  the  attractors  is  highly  enhanced  and  the  regions  in  parameter-space  over  which 
changes  take  place  become  smaller.  In  this  article  we  choose  to  describe  three  additional  windows; 
these  are  the  chaotic  window  entered  immediately  after  the  period-doublings  described  above,  the 
second  is  a  time-periodic  attractor,  and  the  third  Is  a  region  of  parameter  space  that  supports  chaotic, 
steady  and  time-periodic  solutions. 

(vii)  Window  of  Chaotic  Trajectories.  .0252  <  v  $  ,029755. 

No  recognizable  patterns  are  observed  either  in  the  energy  time  series  or  in  the  phase  planes  for 
values  of  the  parameter  in  this  window.  It  should  be  noted,  however,  that  all  the  plotted  quantities 
remain  bounded  in  agreement  with  the  theory  (  Constantin  et  al  [8].)  See  figures  (7.1a-d). 

(viii)  Time  periodic  attractors.  0.023  <>v  <,  Q.0251. 

Another  window  with  time-periodic  attractors  and  at  least  three  period-doubling  bifurcations  was 
also  found.  Unlike  the  previous  periodic  window  the  period  starts  with  small  values  at  the  end  of  the 
window  -  the  period  at  v  =  0.023  is  x  =  .26  -  and  reaches  the  value  x  =  2.86  at  v  =  0.0251  after  three 
period  doublings.  The  results  from  this  window  are  summarized  in  Figures  (8.  la-c)  for  a  representa¬ 
tive  member  v  =  .025. 

(ix)  Window  of  Chaotic  Trajectories.  7  ^  v  £  0.0225. 

Chaotic  solutions  are  also  found  in  this  subwindow  with  no  recognizable  patterns  to  report.  The 
ambiguity  of  the  lower  bound  of  this  window  remains  unresolved  due  to  the  computational  cost  of  the 
required  runs  but  also  due  to  the  well-known  fact  that  this  is  not  the  last  chaotic  window  to  be  found. 
Window  and  subwindow  lengths  become  vanishingly  smaller  thus  halting  any  rational  numerical  pro¬ 
gress.  The  discovery  of  various  attractors  becomes  a  hit  or  miss  problem  and.  for  example,  in  the 
range  .008  SvS  .0225  we  managed  to  find  a  strange  attractor  (v  =  .0225),  a  time-periodic  attractor 
(v  =  .017),  and  bimodal  steady  attractor  (v  =  ,015)  before  another  strange  attractor  is  obtained  at 
v  =  .008.  These  results  are  summarized  in  Figures  (9.1a-e).  For  graphical  clarity,  the  phase-planes  in 
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this  set  of  Figures  were  obtained  by  plotting  (£  .  dEldt )  pairs  which  are  separated  by  a  time  which  is 
a  large  multiple  of  the  time-step. 

S.  Conclusions. 

We  have  carried  out  an  extensive  numerical  sludy  of  the  Kuramoto-Sivashinsky  equation.  Our 
emphasis  is  on  the  characterization  of  transition  to  chaos  and  we  restricted  ourselves  to  odd-parity 
solutions  derived  from  the  same  initial  conditions  and  which  are  easily  reproducible  by  a  single  com¬ 
puter  run.  The  results  for  large  values  of  the  parameter  v,  i.e.  before  chaos  or  multiple  period¬ 
doubling  bifurcations,  are  completely  in  line  with  those  of  previous  investigators  (see  Introduction). 
The  new  behavior  we  found,  and  which  to  our  knowledge  has  not  been  reported  before  for  the 
Kuramoto-Sivashinsky  equation,  is  the  computation  of  several  windows  which  support  time-periodic 
orbits  and  which  themselves  contain  other  subwindows  which  are  separated  by  period-doubling  bifurca¬ 
tions  in  parameter  space.  The  third  periodic  window  we  found  (see  Section  4  above)  provides  substan¬ 
tial  numerical  evidence  of  a  complete  sequence  of  period -doubling  bifurcations  before  chaos  sets  in. 
Based  on  our  numerical  solutions  we  computed  the  Feigenbaum  number  following  the  first  eight  bifur¬ 
cations  of  the  cascade,  and  found  excellent  agreement  with  Feigenbaum's  universal  constant. 
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7.  Tables. 


Overview  of  (lie  attracting  manifolds 

Window 

range 

Description  of  the  attractors 

1 

£ 

V 

< 

OO 

Constant  states. 

.2475 

V 

< 

1 

Fully  modal  steady  attractors. 

.0756 

V 

.2472 

Bimodal  steady  attractors. 

.06697 

£ 

V 

.0755 

Fully  modal  steady  attractors. 

.0599 

£ 

V 

£ 

.06695 

Trimodal  steady  attractors. 

.055238 

V 

<; 

.05985 

Periodic  attractors. 

.03965 

£ 

V 

£ 

.055235 

Fully  modal  steady  attractors. 

.03735 

£ 

V 

£ 

.0396 

Periodic  attractors. 

.0344 

V 

.037348 

Tetramodal  steady  attractors. 

.029756 

V 

£ 

.0343 

Periodic  attractors.  Com¬ 
plete  period-doubling 

sequence  towards  v=.029756. 

.0252 

£ 

V 

.029755 

Chaotic  oscillations. 

.024 

V 

£ 

.0251 

Periodic  attractors.  Com¬ 
plete  period-doubling 

sequence  towards  v=.02512. 

? 

£ 

V 

£ 

0.023 

Chaotic  oscillations. 

Table  1 


First  Periodic  Window  .0598  Svi  .055238 

Viscosity  v 

Period  x 

Figures 

.0598 

1.018 

.0595 

1.03 

.059 

1.05 

.058 

1.08 

.057 

1.13 

(2.1a, b) 

.056 

1.18 

.0559 

1.183 

(2.2a,  b) 

.05585 

2.374 

.0555 

2.48 

(2.3a, b) 

.0553 

2.59 

.05525 

2.663 

.05524 

2.729 

Table  2 
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Third  Periodic  Window  .0343  £  v  <;  .0299756 

Viscosity  v 

Period  x 

Figures 

.0343 

.368 

.0342 

.37 

(6.1a) 

.032 

.44 

.0304 

.89 

.0303 

.89 

(6.1b, c) 

.03007 

.88 

(6.1d) 

.030065 

1.76 

.03006 

1.76 

(6.2a) 

.030015 

1.76 

(6.2b) 

.03 

1.76 

(6.2c) 

.029998 

1.76 

(6.2d) 

.0299975 

3.52 

.029995 

3.52 

(6.3a) 

.02999 

3.52 

(6.3b) 

.029985 

3.52 

.029983 

3.52 

.0299815 

3.52 

(6.3c) 

.029981 

7.04 

.02998 

7.04 

(6.4a) 

.029979 

7.04 

.029978 

7.04 

(6.4b) 

.029976 

7.04 

.0299758 

7.04 

.0299756 

7.04 

(6.4c) 

.029975 

14.08 

(6.5a) 

Tablc  4 
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36-mode  Galerkin  Expansion 

Subwindow  boundary 

Length 

Ratio 

.034615 

mm 

.0303175 

.0042975 

■ 

.0300494 

.0002681 

.02998665 

.00006275 

.029973005 

1.3645xl0_s 

.02997007 

2.935X10-6 

4.65 

.029969435 

6.35x10  7 

4.62 

.029969305 

1.3xl0"7 

4.88 

Table  5 
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8.  Figures 

Figure  (1.1)  :  Steady  attracting  modal  profiles  : 

(a)  v  =  .6  .7  .8  .9  . 

(b)  v~=~.26  .28  .3  .35  .4  .5  . 

Figure  (1.2)  :  Steady  attracting  profiles  : 

(a)  v  =  .2  .22  .24  .248  .28  . 

(b)  v  =  .08  .09  .12  .15  .2  . 

Figure  (1.3)  :  Steady  attracting  trirnodal  profiles  : 

(a)  v  =  .067  .07  .072  .074  .0/55  .08  . 

(b)  v  =  .06  .061  .063  .065  .06695  . 

Figure  (2.1)  :  First  periodic  window,  .0598  £  v  £.055238  : 

(a)  Phase  plane  of  Energy  -  v  =  .057  . 

(b)  Energy  time  series  -  v  =  .057  . 

Figure  (2.2)  :  First  time-periodic  window,  .0598  £  v  £.055238  : 

(a)  Phase  plane  of  Energy  -  v  =  .0559  . 

(b)  Energy  time  series  -  v  =  .0559  . 

Figure  (2.3)  :  First  time-periodic  window,  .0598  £  V  £.055238  : 

(a)  Phase  plane  of  Energy  -  v  =  .0555  . 

(b)  Energy  time  series  -  v  =  .0555  . 

Figure  (3.1)  :  Window  (iii),  fully  modal  steady  attractors,  .06697  £  v  £  .0755 

(a)  Energy  against  time  -  v  =  .0397  . 

(b)  Supremum  of  u(x,t)  against  time  -  v  =  .0397  . 

(c)  Supremum  of  dldx  of  u(jc,r)  against  time  -  v  =  .0397  . 

Figure  (4.1)  :  Second  periodic  window,  time-period  x  =  .68  : 

(a)  Energy  against  time  -  v  =  .038  . 

(b)  Phase  plane  of  Energy  -  v  =  .038  . 

(c)  First  Fourier  coefficient  against  time  -  v  =  .038  . 

(d)  Tenth  Fourier  coefficient  against  time  -  v  =  .038  . 

Figure  (4.2)  :  Second  periodic  window,  time-period  x  =  1.37  : 

(a)  Phase  plane  of  Energy  -  v  =  .03795  . 

(b)  Energy  time  series  -  v  =  .03795  . 

Figure  (4.3)  :  Second  periodic  window,  time-period  x  =  1.37  to  1.43  : 

(a)  Phase  plane  of  Energy  -  v  =  .0379  . 

(b)  Phase  plane  of  Energy  -  v  =  .0377  . 

(c)  Phase  plane  of  the  first  Fourier  coefficient  -  v  =  .0377. 

(d)  Phase  plane  of  the  sixth  Fourier  coefficient  -  v  =  .0377. 
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Figure  (4.4)  :  Second  periodic  window,  time-period  x  =  2.87  : 

(a)  Phase  plane  of  Energy  -  v  =  .0376  . 

(b)  Time  series  of  sixth  Fourier  coefficient  -  v  =  .0376  . 

(c)  Phase  plane  of  die  sixth  Fourier  coefficient  -  v  =  .0376. 

Figure  (4.5)  :  Second  periodic  window,  period-halving  from  x  =  3.01  to  1.51 

(a)  Energy  time  series  -  v  =  .03736  . 

Phase  plane  of  Energy  -  v  =  .03736  . 

Phase  plane  of  Energy  -  v  =  .03735  . 

Figure  (4.6)  :  Second  periodic  window  : 

Evolution  of  the  spatial  profile  in  the  periodic  attractor  -  v  =  .039. 

Figure  (4.7)  :  Second  periodic  window  : 

Evolution  of  die  spatial  profile  in  the  periodic  attractor  -  v  =  .0375. 

Figure  (5.1)  :  Window  (v),  steady  tetramodal  attractors,  .0344  ^  v  <  .037348 

(a)  Energy  time  series  -  v  =  .0373  . 

(b)  Steady  tetramodal  attracting  profile  -  v  =  .0373  . 

Figure  (6.1)  :  Third  periodic  window,  first  period -doubling  : 

(a)  Phase  plane  of  Energy  -  v  =  .0342  . 

(b)  Phase  plane  of  Energy  -  v  =  .0303  . 

(c)  Energy  time  series  -  v  =  .0303  . 

(d)  Phase  plane  of  Energy  -  v  =  .03007  . 

Figure  (6.2)  :  Third  periodic  window,  second  period-doubling  : 

(a)  Phase  plane  of  Energy  -  v  =  .03006  . 

(b)  Phase  plane  of  Energy  -  v  =  .030015  . 

(c)  Phase  plane  of  Energy  -  v  =  .03  . 

(d)  Phase  plane  of  Energy  -  v  =  .029998  . 

Figure  (6.3)  :  Third  periodic  window,  third  period-doubling  : 

(a)  Phase  plane  of  Energy  -  v  =  .029995  . 

(b)  Phase  plane  of  Energy  -  v  =  .02999  . 

(c)  Phase  plane  of  Energy  -  v  =  .0299815  . 

Figure  (6.4)  :  Third  periodic  window,  fourth  period-doubling  : 

(a)  Phase  plane  of  Energy  -  v  =  .02998  . 

(b)  Phase  plane  of  Energy  -  v  =  .029978  . 

(c)  Phase  plane  of  Energy  -  v  =  .0299756  . 

Figure  (6.5)  :  Third  periodic  window,  fifth  period-doubling  : 

(a)  Phase  plane  of  Energy  -  v  =  .029975  . 

(b)  Phase  plane  of  Energy  -  v  =  .02997  . 

(c)  Phase  plane  of  Energy  -  v  =  .02995  . 
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(d)  Phase  plane  of  Energy  -  v  =  .0299  . 

Figure  (7.1)  :  Chaotic  windows,  .0252  .029755  : 

(a)  Energy  time  series  (Chaotic  case)  -  v  =  .028  . 

(b)  Energy  time  series  (Chaotic  case)  -  v  =  .027  . 

(c)  Energy  time  series  (Chaotic  case)  -  v  =  .0255  . 

(d)  Energy  time  series  (Chaotic  case)  -  v  =  .0252  . 

Figure  (8.1)  :  Time-periodic  attractors,  .023  .0251  : 

(a)  Energy  time  series  (Periodic  case)  -  v  =  .025  . 

(b)  Energy  time  series  (Blow-up)  -  v  =  .025  . 

(c)  Phase  plane  of  Energy  -  v  =  .025  . 

Figure  (9.1)  :  Chaotic  windows,  ?  <,  v  <,  .0255  : 

(a)  Phase  plane  of  Energy  (Chaotic  case)  -  v  =  .0225  . 

(b)  Phase  plane  of  Energy  (Periodic  case)  -  v  =  .017  . 

(c)  Energy  time  series  (Periodic  case)  -  v  =  .017  . 

(d)  Energy  time  series  (Bimodal  steady  attractor)  -  v  =  .015  . 

(e)  Phase  plane  of  Energy  (Chaotic  case)  -  v  =  .008  . 
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